#ifndef __POLYNOMICON__
#define __POLYNOMICON__


#include "cxfloat.h"
#include <vector>
#include "Pincushion.h"


class FourierTron;


class Polynomicon
{
  public:
    Polynomicon(FourierTron* f, std::vector<SeamStress::Seamstress*>* ss, unsigned long int nt=1, mpfr_prec_t p=4);
    ~Polynomicon();
    
    void calculateFactorials(unsigned long int N);
    void factorialVal(unsigned long int ind, refloat& val);
    void factorialRecipVal(unsigned long int ind, refloat& val);
    void factorialVal(unsigned long int ind, cxfloat& val);
    void factorialRecipVal(unsigned long int ind, cxfloat& val);
    
    void multiplyPoly(const std::vector<cxfloat>& in1, const std::vector<cxfloat>& in2, std::vector<cxfloat>& out);
    void multiplyPoly(const refloat_vec& in1, const refloat_vec& in2, refloat_vec& out);
    void multiplyPoly(const refloat_vec& in1, const refloat_vec& in2, refloat_vec& out, long int log2scale);
    
    void shiftPoly(std::vector<cxfloat>& poly, const cxfloat& Delta, bool set_scale=false, long int log_scale=0);
    void shiftPolyBrute(std::vector<cxfloat>& poly, const cxfloat& Delta);
    void shiftPoly(refloat_vec& poly, const refloat& Delta, bool set_scale=false, long int log_scale=0);
    void shiftPolyBrute(refloat_vec& poly, const refloat& Delta);
    
    void scalePoly(std::vector<cxfloat>& poly, const cxfloat& scale);
    void scalePoly(refloat_vec& poly, const refloat& scale);
    
    void vectorMultiply(std::vector<cxfloat>& poly, const cxfloat& scale);
    void vectorMultiply(refloat_vec& poly, const refloat& scale);
    
    void reciprocalPowerSeries(unsigned long int nterms, const std::vector<cxfloat>& poly, std::vector<cxfloat>& series);
    void reciprocalPowerSeries(unsigned long int nterms, const refloat_vec& poly, refloat_vec& series);
    
    void addPoly(const std::vector<cxfloat>& to_add, std::vector<cxfloat>& added_to);
    void addPoly(const refloat_vec& to_add, refloat_vec& added_to);
    void subPoly(const std::vector<cxfloat>& to_subtract, std::vector<cxfloat>& subtracted_from);
    
    void evaluatePoly(refloat& output, const refloat_vec& poly, const refloat& input);
    
    void generalizedShift(refloat_vec& vec1, refloat_vec& vec2, refloat_vec& vec3, refloat_vec& output, long int log2scale);
    
    void composePolyRecursive(unsigned long int n, const refloat_vec& Q, const refloat_vec& P, refloat_vec& R);
    void composePoly(const refloat_vec& Q, const refloat_vec& P, refloat_vec& R);
    void differentiatePoly(const refloat_vec& input, refloat_vec& output);
    
    void composeQuadraticPoly(const refloat_vec& Q, const refloat& a, const refloat& b, refloat_vec& R);
    
    void modPoly(const refloat_vec& input, const refloat_vec& divisor, refloat_vec& remainder);
    void modPoly(const refloat_vec& input, const refloat_vec& divisor, const refloat_vec& rev_b_inv, refloat_vec& remainder);
    
    
    FourierTron* getFourierTron();
    
  protected:
    SeamStress::Pincushion<Polynomicon> pins;
    
    //floating point representation of factorial values
    refloat_vec factorial;
    //floating point representation of reciprocal factorial values
    refloat_vec factorial_inv;
    
    FourierTron* fourier;
    
    mp_prec_t prec;
    
    unsigned long int nthreads;
    
  private:
    
    void factorialInit(void* arg);
    void factorialRun(void* arg);
    void factorialInv(void* arg);
    
    void multiplyPolyThread1(void* arg);
    void multiplyPolySchoolbook(void* arg);
    
    void scalePolyPower2(void* arg);
    void setUPoly(void* arg);
    void setVPoly(void* arg);
    void finalizeShift(void* arg);
    
    void scalePolyPower2Real(void* arg);
    void setUPolyReal(void* arg);
    void setVPolyReal(void* arg);
    void finalizeShiftReal(void* arg);
    
    void powerUp(void* arg);
    void powerUpReal(void* arg);
    
    void vectorMultiply(void* arg);
    void vectorMultiplyReal(void* arg);
    
    void addPolyThread1(void* arg);
    void addPolyThread1Real(void* arg);
    void subPolyThread1(void* arg);
    
    void addIntermediate(void* arg);
    void parallelHornerReal(void* arg);
    
    void reflectAboutCenter(void* arg);
    void oneToOneMultiply(void* arg);
    
    void differentiatePolyThreadReal(void* arg);
    
    void modPolyThread1(void* arg);
    
    unsigned long int thread_tot;
    unsigned long int thread_var1;
    unsigned long int thread_var2;
    unsigned long int thread_var3;
    std::vector<cxfloat> thread_vec1;//one entry per thread
    refloat_vec rethread_vec1;//one entry per thread
    refloat_vec rethread_vec2;
    
    cxfloat tempval1, tempval2;
    refloat retempval1, retempval2;
    
    std::vector<cxfloat> tempvec1;
    std::vector<cxfloat> tempvec2;
    std::vector<cxfloat> tempvec3;
    
    refloat_vec retempvec1;
    refloat_vec retempvec2;
    refloat_vec retempvec3;
    
    std::vector<cxfloat>* tempvecpointer1;
    const std::vector<cxfloat>* tempvecpointer2;
    refloat_vec* retempvecpointer1;
    const refloat_vec* retempvecpointer2;
    const refloat_vec* retempvecpointer3;
    
    long int p2_scale;
};


inline void Polynomicon::factorialVal(unsigned long int ind, refloat& val){val.set(factorial[ind]);}
inline void Polynomicon::factorialRecipVal(unsigned long int ind, refloat& val){val.set(factorial_inv[ind]);}
inline void Polynomicon::factorialVal(unsigned long int ind, cxfloat& val){retempval1.set(0);val.set(factorial[ind], retempval1);}
inline void Polynomicon::factorialRecipVal(unsigned long int ind, cxfloat& val){retempval1.set(0);val.set(factorial_inv[ind], retempval1);}

inline FourierTron* Polynomicon::getFourierTron(){return fourier;}


#endif
